1st part

ChatGPT prompt: How to plot the histogram of the pixel value distribution?

from matplotlib import pyplot as plt
import numpy as np
from PIL import Image
import seaborn as sns
from scipy.stats import norm

# Load the image
image_path = './img.jpg'
image = Image.open(image_path)

# Convert image to grayscale
gray_image = image.convert("L")

# Convert grayscale image to a numpy array
image_array = np.array(gray_image)

# Flatten the array to 1D
pixel_values = image_array.flatten()

# Plot histogram of pixel values
plt.figure(figsize=(10, 6))
sns.histplot(pixel_values, bins=256, kde=False, color='blue')

# Fit a normal distribution to the pixel values
mean, std_dev = norm.fit(pixel_values)

# Plot the PDF of the normal distribution using the mean and standard deviation obtained above
xmin, xmax = plt.xlim()
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, mean, std_dev)

plt.plot(x, p, 'k', linewidth=2)
title = "Fit results: mean = %.2f,  std = %.2f" % (mean, std_dev)
plt.title(title)

plt.xlabel('Pixel Intensity')
plt.ylabel('Frequency')

# Display the histogram
plt.show()

# Return mean and standard deviation for further analysis
mean, std_dev

(82.35787963867188, 9.168225715863297)

ChatGPT prompt: How to plot the control chart of this modified image?

# Calculate the bounds for the lower and upper 0.001 probability limits (quantiles)
lower_bound = norm.ppf(0.001, mean, std_dev)
upper_bound = norm.ppf(0.999, mean, std_dev)

# Identify the pixels that are out of these bounds
outliers = np.where((pixel_values < lower_bound) | (pixel_values > upper_bound))

# Create a copy of the image array for modification
modified_image_array = image_array.copy()

# Set the outlier pixel values to zero (black)
modified_image_array.flat[outliers] = 0

# Convert the modified array back to an image
modified_image = Image.fromarray(modified_image_array.astype(np.uint8))

# Plotting the original and the modified image side by side
fig, axs = plt.subplots(2, 2, figsize=(12, 10))

# Original image
axs[0, 0].imshow(image_array, cmap='gray')
axs[0, 0].set_title('Original Image')
axs[0, 0].axis('off')

# Modified image with outliers set to black
axs[0, 1].imshow(modified_image_array, cmap='gray')
axs[0, 1].set_title('Modified Image with Outliers')
axs[0, 1].axis('off')

# Control chart
# Normal points
axs[1, 0].plot(pixel_values, marker='o', linestyle='', markersize=0.5, color='blue')
# Out of control points
axs[1, 0].plot(outliers[0], pixel_values[outliers], marker='o', linestyle='', markersize=0.5, color='red')
# Mean, lower and upper bounds
axs[1, 0].axhline(y=mean, color='r', linestyle='-', label='Mean')
axs[1, 0].axhline(y=lower_bound, color='g', linestyle='--', label='Lower Bound')
axs[1, 0].axhline(y=upper_bound, color='b', linestyle='--', label='Upper Bound')
axs[1, 0].set_title('Control Chart of Pixel Values')
axs[1, 0].legend()

# Hide the fourth subplot (if not used)
axs[1, 1].axis('off')

# Display the plots
plt.tight_layout()
plt.show()

# Outputting the bounds and the number of identified outliers
(lower_bound, upper_bound, len(outliers[0]))

(54.02593234127259, 110.68982693607116, 846)

len(outliers[0])

846

The ChatGPT prompt of this part was a little bit complicated, we tried to teach the process patch approach first and then ask for the following code.

# Define the window size
window_size = 51

# Function to process each patch
def process_patch(patch, mean, std_dev):
    # Flatten the patch to 1D array of pixel values
    patch_values = patch.flatten()
    # Calculate the bounds for the 0.001 probability limits
    lower_bound = norm.ppf(0.001, mean, std_dev)
    upper_bound = norm.ppf(0.999, mean, std_dev)
    # Identify pixels that are out of these bounds
    outliers = np.where((patch_values < lower_bound) | (patch_values > upper_bound))
    # Return the indexes of the outliers within the patch
    return outliers[0]

# Prepare an empty array to store the modified image
modified_image_patches = np.copy(image_array)

# Loop over the image array by the defined window size
for i in range(0, image_array.shape[0] - window_size + 1, window_size):
    for j in range(0, image_array.shape[1] - window_size + 1, window_size):
        # Extract the patch
        patch = image_array[i:i+window_size, j:j+window_size]
        # Fit a normal distribution to the pixel values of the patch
        patch_mean, patch_std_dev = norm.fit(patch.flatten())
        # Identify outliers in the patch
        outlier_indexes = process_patch(patch, patch_mean, patch_std_dev)
        # If outliers are found, set their values to zero (black) in the corresponding positions
        if len(outlier_indexes) > 0:
            # Convert the 1D outlier indexes back to 2D indexes
            outlier_indexes_2d = np.unravel_index(outlier_indexes, (window_size, window_size))
            # Set the outliers to zero in the modified image
            modified_image_patches[i+outlier_indexes_2d[0], j+outlier_indexes_2d[1]] = 0

# Convert the modified array back to an image
modified_image_patches = Image.fromarray(modified_image_patches)

# Plot the original and the modified image side by side
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# Original image
ax1.imshow(image_array, cmap='gray')
ax1.set_title('Original Image')
ax1.axis('off')

# Modified image with outliers set to black based on patches
ax2.imshow(modified_image_patches, cmap='gray')
ax2.set_title('Modified Image with Patch Outliers')
ax2.axis('off')

# Display the plots
plt.tight_layout()
plt.show()

len(outliers[0])

846

2nd part

ChatGPT prompt: For each row and column (seperately) of your image, construct an appropriate chart for monitoring the mean and variance. Find the pixels that are out of control. After finding those pixels, change the value of these pixels to zero.

# Re-importing necessary libraries as the execution state has been reset.
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt

# Re-load the image
image = Image.open("img.jpg")
image_array = np.array(gray_image)

# Function to identify out-of-control pixels based on mean for rows or columns
def identify_out_of_control_on_mean(data):
    # Calculate the mean and standard deviation for each row or column
    mean = np.mean(data, axis=1)
    std_dev = np.std(data, axis=1)

    # Calculate the control limits for mean (assuming a normal distribution)
    # We will use the 3-sigma rule for control limits
    upper_control_limit_mean = mean + 3 * std_dev
    lower_control_limit_mean = mean - 3 * std_dev

    # Identify the out-of-control pixels for mean
    out_of_control_mean = np.where((data > upper_control_limit_mean[:, None]) |
                                   (data < lower_control_limit_mean[:, None]))

    return out_of_control_mean

# Identify out-of-control pixels for rows based on mean
out_of_control_rows = identify_out_of_control_on_mean(image_array)

# Create a copy of the image array to modify pixels that are out of control for rows
modified_image_rows_mean = np.copy(image_array)
modified_image_rows_mean[out_of_control_rows] = 0

# Identify out-of-control pixels for columns based on mean
out_of_control_columns = identify_out_of_control_on_mean(image_array.T)

# Create another copy of the image array to modify pixels that are out of control for columns
modified_image_columns_mean = np.copy(image_array)
modified_image_columns_mean[:, out_of_control_columns[0]] = 0

# Convert the modified arrays back to images
modified_image_rows_mean = Image.fromarray(modified_image_rows_mean)
modified_image_columns_mean = Image.fromarray(modified_image_columns_mean)

# Plot the original and the modified images for rows and columns based on mean
fig, axs = plt.subplots(2, 2, figsize=(12, 12))

# Original image
axs[0, 0].imshow(image_array, cmap='gray')
axs[0, 0].set_title('Original Image')
axs[0, 0].axis('off')

# Modified image with rows processed based on mean
axs[0, 1].imshow(modified_image_rows_mean, cmap='gray')
axs[0, 1].set_title('Rows Processed on Mean')
axs[0, 1].axis('off')

# Original image for symmetry in layout
axs[1, 0].imshow(image_array, cmap='gray')
axs[1, 0].set_title('Original Image')
axs[1, 0].axis('off')

# Modified image with columns processed based on mean
axs[1, 1].imshow(modified_image_columns_mean, cmap='gray')
axs[1, 1].set_title('Columns Processed on Mean')
axs[1, 1].axis('off')

# Display the plots
plt.tight_layout()
plt.show()

len(out_of_control_rows[0])+len(out_of_control_rows[1])

1776

len(out_of_control_columns[0])+len(out_of_control_columns[1])

1882

Our Proposal Moving Average Filtering

import numpy as np
from scipy.ndimage import uniform_filter

# Define the moving average filter function
def moving_average_filter(image, kernel_size=3):
    """
    Apply a moving average filter to an image.

    Parameters:
    - image: A 2D numpy array representing the image.
    - kernel_size: The size of the kernel (window) for the moving average. Default is 3 .

    Returns:
    - A 2D numpy array representing the filtered image.
    """
    # Apply the uniform filter, which is a moving average filter
    filtered_image = uniform_filter(image, size=kernel_size)
    return filtered_image

# Load the image and convert it to grayscale
image_path = './img.jpg'
image = Image.open(image_path).convert('L')
image_array = np.array(gray_image)

# Apply the moving average filter with a kernel size of 5
filtered_image_array = moving_average_filter(image_array, kernel_size=5)

# Convert the filtered array back to an image for display
filtered_image = Image.fromarray(filtered_image_array.astype(np.uint8))

# Flatten the original and filtered image arrays for analysis
original_pixel_values = image_array.flatten()
filtered_pixel_values = filtered_image_array.flatten()

# Calculate mean and standard deviation for both original and filtered images
original_mean = np.mean(original_pixel_values)
original_std_dev = np.std(original_pixel_values)
filtered_mean = np.mean(filtered_pixel_values)
filtered_std_dev = np.std(filtered_pixel_values)

# Calculate control limits for both images
original_lower_bound = norm.ppf(0.001, original_mean, original_std_dev)
original_upper_bound = norm.ppf(0.999, original_mean, original_std_dev)
filtered_lower_bound = norm.ppf(0.001, filtered_mean, filtered_std_dev)
filtered_upper_bound = norm.ppf(0.999, filtered_mean, filtered_std_dev)

# Identify outliers for both images
original_outliers = np.where((original_pixel_values < original_lower_bound) | (original_pixel_values > original_upper_bound))
filtered_outliers = np.where((filtered_pixel_values < filtered_lower_bound) | (filtered_pixel_values > filtered_upper_bound))

# Plotting
fig, axs = plt.subplots(2, 2, figsize=(15, 12))

# Original grayscale image and its control chart
axs[0, 0].imshow(image_array, cmap='gray')
axs[0, 0].set_title('Original Grayscale Image')
axs[0, 0].axis('off')
axs[1, 0].plot(original_pixel_values, marker='o', linestyle='', markersize=0.5, color='blue')
axs[1, 0].plot(original_outliers[0], original_pixel_values[original_outliers], marker='o', linestyle='', markersize=0.5, color='red')
axs[1, 0].axhline(y=original_mean, color='r', linestyle='-', label='Mean')
axs[1, 0].axhline(y=original_lower_bound, color='g', linestyle='--', label='Lower Bound')
axs[1, 0].axhline(y=original_upper_bound, color='b', linestyle='--', label='Upper Bound')
axs[1, 0].set_title('Control Chart for Original Image')
axs[1, 0].legend()

# Filtered image and its control chart
axs[0, 1].imshow(filtered_image, cmap='gray')
axs[0, 1].set_title('Moving Average Filtered Image')
axs[0, 1].axis('off')
axs[1, 1].plot(filtered_pixel_values, marker='o', linestyle='', markersize=0.5, color='blue')
axs[1, 1].plot(filtered_outliers[0], filtered_pixel_values[filtered_outliers], marker='o', linestyle='', markersize=0.5, color='red')
axs[1, 1].axhline(y=filtered_mean, color='r', linestyle='-', label='Mean')
axs[1, 1].axhline(y=filtered_lower_bound, color='g', linestyle='--', label='Lower Bound')
axs[1, 1].axhline(y=filtered_upper_bound, color='b', linestyle='--', label='Upper Bound')
axs[1, 1].set_title('Control Chart for Filtered Image')
axs[1, 1].legend()

# Display the plots
plt.tight_layout()
plt.show()

len(original_outliers[0])

846

len(filtered_outliers[0])

563

Alternative images

import random as rd
rd.seed(177)

rd.randint(2,196)

#109
#151
#183
#105
#13 images are choosen

13

import numpy as np
from scipy.ndimage import uniform_filter

# Define the moving average filter function
def moving_average_filter(image, kernel_size=3):
    """
    Apply a moving average filter to an image.

    Parameters:
    - image: A 2D numpy array representing the image.
    - kernel_size: The size of the kernel (window) for the moving average. Default is 3 .

    Returns:
    - A 2D numpy array representing the filtered image.
    """
    # Apply the uniform filter, which is a moving average filter
    filtered_image = uniform_filter(image, size=kernel_size)
    return filtered_image

# Load the image and convert it to grayscale
image_path = './0013.jpg'
image = Image.open(image_path).convert('L')
image_array = np.array(image)

# Apply the moving average filter with a kernel size of 5
filtered_image_array = moving_average_filter(image_array, kernel_size=5)

# Convert the filtered array back to an image for display
filtered_image = Image.fromarray(filtered_image_array.astype(np.uint8))

# Flatten the original and filtered image arrays for analysis
original_pixel_values = image_array.flatten()
filtered_pixel_values = filtered_image_array.flatten()

# Calculate mean and standard deviation for both original and filtered images
original_mean = np.mean(original_pixel_values)
original_std_dev = np.std(original_pixel_values)
filtered_mean = np.mean(filtered_pixel_values)
filtered_std_dev = np.std(filtered_pixel_values)

# Calculate control limits for both images
original_lower_bound = norm.ppf(0.001, original_mean, original_std_dev)
original_upper_bound = norm.ppf(0.999, original_mean, original_std_dev)
filtered_lower_bound = norm.ppf(0.001, filtered_mean, filtered_std_dev)
filtered_upper_bound = norm.ppf(0.999, filtered_mean, filtered_std_dev)

# Identify outliers for both images
original_outliers = np.where((original_pixel_values < original_lower_bound) | (original_pixel_values > original_upper_bound))
filtered_outliers = np.where((filtered_pixel_values < filtered_lower_bound) | (filtered_pixel_values > filtered_upper_bound))

# Plotting
fig, axs = plt.subplots(2, 2, figsize=(15, 12))

# Original grayscale image and its control chart
axs[0, 0].imshow(image_array, cmap='gray')
axs[0, 0].set_title('Original Grayscale Image')
axs[0, 0].axis('off')
axs[1, 0].plot(original_pixel_values, marker='o', linestyle='', markersize=0.5, color='blue')
axs[1, 0].plot(original_outliers[0], original_pixel_values[original_outliers], marker='o', linestyle='', markersize=0.5, color='red')
axs[1, 0].axhline(y=original_mean, color='r', linestyle='-', label='Mean')
axs[1, 0].axhline(y=original_lower_bound, color='g', linestyle='--', label='Lower Bound')
axs[1, 0].axhline(y=original_upper_bound, color='b', linestyle='--', label='Upper Bound')
axs[1, 0].set_title('Control Chart for Original Image')
axs[1, 0].legend()

# Filtered image and its control chart
axs[0, 1].imshow(filtered_image, cmap='gray')
axs[0, 1].set_title('Moving Average Filtered Image')
axs[0, 1].axis('off')
axs[1, 1].plot(filtered_pixel_values, marker='o', linestyle='', markersize=0.5, color='blue')
axs[1, 1].plot(filtered_outliers[0], filtered_pixel_values[filtered_outliers], marker='o', linestyle='', markersize=0.5, color='red')
axs[1, 1].axhline(y=filtered_mean, color='r', linestyle='-', label='Mean')
axs[1, 1].axhline(y=filtered_lower_bound, color='g', linestyle='--', label='Lower Bound')
axs[1, 1].axhline(y=filtered_upper_bound, color='b', linestyle='--', label='Upper Bound')
axs[1, 1].set_title('Control Chart for Filtered Image')
axs[1, 1].legend()

# Display the plots
plt.tight_layout()
plt.show()

len(original_outliers[0]), len(filtered_outliers[0])

(50, 19)

import numpy as np
from scipy.ndimage import uniform_filter

# Define the moving average filter function
def moving_average_filter(image, kernel_size=3):
    """
    Apply a moving average filter to an image.

    Parameters:
    - image: A 2D numpy array representing the image.
    - kernel_size: The size of the kernel (window) for the moving average. Default is 3 .

    Returns:
    - A 2D numpy array representing the filtered image.
    """
    # Apply the uniform filter, which is a moving average filter
    filtered_image = uniform_filter(image, size=kernel_size)
    return filtered_image

# Load the image and convert it to grayscale
image_path = './0105.jpg'
image = Image.open(image_path).convert('L')
image_array = np.array(image)

# Apply the moving average filter with a kernel size of 5
filtered_image_array = moving_average_filter(image_array, kernel_size=5)

# Convert the filtered array back to an image for display
filtered_image = Image.fromarray(filtered_image_array.astype(np.uint8))

# Flatten the original and filtered image arrays for analysis
original_pixel_values = image_array.flatten()
filtered_pixel_values = filtered_image_array.flatten()

# Calculate mean and standard deviation for both original and filtered images
original_mean = np.mean(original_pixel_values)
original_std_dev = np.std(original_pixel_values)
filtered_mean = np.mean(filtered_pixel_values)
filtered_std_dev = np.std(filtered_pixel_values)

# Calculate control limits for both images
original_lower_bound = norm.ppf(0.001, original_mean, original_std_dev)
original_upper_bound = norm.ppf(0.999, original_mean, original_std_dev)
filtered_lower_bound = norm.ppf(0.001, filtered_mean, filtered_std_dev)
filtered_upper_bound = norm.ppf(0.999, filtered_mean, filtered_std_dev)

# Identify outliers for both images
original_outliers = np.where((original_pixel_values < original_lower_bound) | (original_pixel_values > original_upper_bound))
filtered_outliers = np.where((filtered_pixel_values < filtered_lower_bound) | (filtered_pixel_values > filtered_upper_bound))

# Plotting
fig, axs = plt.subplots(2, 2, figsize=(15, 12))

# Original grayscale image and its control chart
axs[0, 0].imshow(image_array, cmap='gray')
axs[0, 0].set_title('Original Grayscale Image')
axs[0, 0].axis('off')
axs[1, 0].plot(original_pixel_values, marker='o', linestyle='', markersize=0.5, color='blue')
axs[1, 0].plot(original_outliers[0], original_pixel_values[original_outliers], marker='o', linestyle='', markersize=0.5, color='red')
axs[1, 0].axhline(y=original_mean, color='r', linestyle='-', label='Mean')
axs[1, 0].axhline(y=original_lower_bound, color='g', linestyle='--', label='Lower Bound')
axs[1, 0].axhline(y=original_upper_bound, color='b', linestyle='--', label='Upper Bound')
axs[1, 0].set_title('Control Chart for Original Image')
axs[1, 0].legend()

# Filtered image and its control chart
axs[0, 1].imshow(filtered_image, cmap='gray')
axs[0, 1].set_title('Moving Average Filtered Image')
axs[0, 1].axis('off')
axs[1, 1].plot(filtered_pixel_values, marker='o', linestyle='', markersize=0.5, color='blue')
axs[1, 1].plot(filtered_outliers[0], filtered_pixel_values[filtered_outliers], marker='o', linestyle='', markersize=0.5, color='red')
axs[1, 1].axhline(y=filtered_mean, color='r', linestyle='-', label='Mean')
axs[1, 1].axhline(y=filtered_lower_bound, color='g', linestyle='--', label='Lower Bound')
axs[1, 1].axhline(y=filtered_upper_bound, color='b', linestyle='--', label='Upper Bound')
axs[1, 1].set_title('Control Chart for Filtered Image')
axs[1, 1].legend()

# Display the plots
plt.tight_layout()
plt.show()

len(original_outliers[0]), len(filtered_outliers[0])

(164, 1107)

import numpy as np
from scipy.ndimage import uniform_filter

# Define the moving average filter function
def moving_average_filter(image, kernel_size=3):
    """
    Apply a moving average filter to an image.

    Parameters:
    - image: A 2D numpy array representing the image.
    - kernel_size: The size of the kernel (window) for the moving average. Default is 3 .

    Returns:
    - A 2D numpy array representing the filtered image.
    """
    # Apply the uniform filter, which is a moving average filter
    filtered_image = uniform_filter(image, size=kernel_size)
    return filtered_image

# Load the image and convert it to grayscale
image_path = './0109.jpg'
image = Image.open(image_path).convert('L')
image_array = np.array(image)

# Apply the moving average filter with a kernel size of 5
filtered_image_array = moving_average_filter(image_array, kernel_size=5)

# Convert the filtered array back to an image for display
filtered_image = Image.fromarray(filtered_image_array.astype(np.uint8))

# Flatten the original and filtered image arrays for analysis
original_pixel_values = image_array.flatten()
filtered_pixel_values = filtered_image_array.flatten()

# Calculate mean and standard deviation for both original and filtered images
original_mean = np.mean(original_pixel_values)
original_std_dev = np.std(original_pixel_values)
filtered_mean = np.mean(filtered_pixel_values)
filtered_std_dev = np.std(filtered_pixel_values)

# Calculate control limits for both images
original_lower_bound = norm.ppf(0.001, original_mean, original_std_dev)
original_upper_bound = norm.ppf(0.999, original_mean, original_std_dev)
filtered_lower_bound = norm.ppf(0.001, filtered_mean, filtered_std_dev)
filtered_upper_bound = norm.ppf(0.999, filtered_mean, filtered_std_dev)

# Identify outliers for both images
original_outliers = np.where((original_pixel_values < original_lower_bound) | (original_pixel_values > original_upper_bound))
filtered_outliers = np.where((filtered_pixel_values < filtered_lower_bound) | (filtered_pixel_values > filtered_upper_bound))

# Plotting
fig, axs = plt.subplots(2, 2, figsize=(15, 12))

# Original grayscale image and its control chart
axs[0, 0].imshow(image_array, cmap='gray')
axs[0, 0].set_title('Original Grayscale Image')
axs[0, 0].axis('off')
axs[1, 0].plot(original_pixel_values, marker='o', linestyle='', markersize=0.5, color='blue')
axs[1, 0].plot(original_outliers[0], original_pixel_values[original_outliers], marker='o', linestyle='', markersize=0.5, color='red')
axs[1, 0].axhline(y=original_mean, color='r', linestyle='-', label='Mean')
axs[1, 0].axhline(y=original_lower_bound, color='g', linestyle='--', label='Lower Bound')
axs[1, 0].axhline(y=original_upper_bound, color='b', linestyle='--', label='Upper Bound')
axs[1, 0].set_title('Control Chart for Original Image')
axs[1, 0].legend()

# Filtered image and its control chart
axs[0, 1].imshow(filtered_image, cmap='gray')
axs[0, 1].set_title('Moving Average Filtered Image')
axs[0, 1].axis('off')
axs[1, 1].plot(filtered_pixel_values, marker='o', linestyle='', markersize=0.5, color='blue')
axs[1, 1].plot(filtered_outliers[0], filtered_pixel_values[filtered_outliers], marker='o', linestyle='', markersize=0.5, color='red')
axs[1, 1].axhline(y=filtered_mean, color='r', linestyle='-', label='Mean')
axs[1, 1].axhline(y=filtered_lower_bound, color='g', linestyle='--', label='Lower Bound')
axs[1, 1].axhline(y=filtered_upper_bound, color='b', linestyle='--', label='Upper Bound')
axs[1, 1].set_title('Control Chart for Filtered Image')
axs[1, 1].legend()

# Display the plots
plt.tight_layout()
plt.show()

len(original_outliers[0]), len(filtered_outliers[0])

(1503, 720)

import numpy as np
from scipy.ndimage import uniform_filter

# Define the moving average filter function
def moving_average_filter(image, kernel_size=3):
    """
    Apply a moving average filter to an image.

    Parameters:
    - image: A 2D numpy array representing the image.
    - kernel_size: The size of the kernel (window) for the moving average. Default is 3 .

    Returns:
    - A 2D numpy array representing the filtered image.
    """
    # Apply the uniform filter, which is a moving average filter
    filtered_image = uniform_filter(image, size=kernel_size)
    return filtered_image

# Load the image and convert it to grayscale
image_path = './0151.jpg'
image = Image.open(image_path).convert('L')
image_array = np.array(image)

# Apply the moving average filter with a kernel size of 5
filtered_image_array = moving_average_filter(image_array, kernel_size=5)

# Convert the filtered array back to an image for display
filtered_image = Image.fromarray(filtered_image_array.astype(np.uint8))

# Flatten the original and filtered image arrays for analysis
original_pixel_values = image_array.flatten()
filtered_pixel_values = filtered_image_array.flatten()

# Calculate mean and standard deviation for both original and filtered images
original_mean = np.mean(original_pixel_values)
original_std_dev = np.std(original_pixel_values)
filtered_mean = np.mean(filtered_pixel_values)
filtered_std_dev = np.std(filtered_pixel_values)

# Calculate control limits for both images
original_lower_bound = norm.ppf(0.001, original_mean, original_std_dev)
original_upper_bound = norm.ppf(0.999, original_mean, original_std_dev)
filtered_lower_bound = norm.ppf(0.001, filtered_mean, filtered_std_dev)
filtered_upper_bound = norm.ppf(0.999, filtered_mean, filtered_std_dev)

# Identify outliers for both images
original_outliers = np.where((original_pixel_values < original_lower_bound) | (original_pixel_values > original_upper_bound))
filtered_outliers = np.where((filtered_pixel_values < filtered_lower_bound) | (filtered_pixel_values > filtered_upper_bound))

# Plotting
fig, axs = plt.subplots(2, 2, figsize=(15, 12))

# Original grayscale image and its control chart
axs[0, 0].imshow(image_array, cmap='gray')
axs[0, 0].set_title('Original Grayscale Image')
axs[0, 0].axis('off')
axs[1, 0].plot(original_pixel_values, marker='o', linestyle='', markersize=0.5, color='blue')
axs[1, 0].plot(original_outliers[0], original_pixel_values[original_outliers], marker='o', linestyle='', markersize=0.5, color='red')
axs[1, 0].axhline(y=original_mean, color='r', linestyle='-', label='Mean')
axs[1, 0].axhline(y=original_lower_bound, color='g', linestyle='--', label='Lower Bound')
axs[1, 0].axhline(y=original_upper_bound, color='b', linestyle='--', label='Upper Bound')
axs[1, 0].set_title('Control Chart for Original Image')
axs[1, 0].legend()

# Filtered image and its control chart
axs[0, 1].imshow(filtered_image, cmap='gray')
axs[0, 1].set_title('Moving Average Filtered Image')
axs[0, 1].axis('off')
axs[1, 1].plot(filtered_pixel_values, marker='o', linestyle='', markersize=0.5, color='blue')
axs[1, 1].plot(filtered_outliers[0], filtered_pixel_values[filtered_outliers], marker='o', linestyle='', markersize=0.5, color='red')
axs[1, 1].axhline(y=filtered_mean, color='r', linestyle='-', label='Mean')
axs[1, 1].axhline(y=filtered_lower_bound, color='g', linestyle='--', label='Lower Bound')
axs[1, 1].axhline(y=filtered_upper_bound, color='b', linestyle='--', label='Upper Bound')
axs[1, 1].set_title('Control Chart for Filtered Image')
axs[1, 1].legend()

# Display the plots
plt.tight_layout()
plt.show()

len(original_outliers[0]), len(filtered_outliers[0])

(832, 813)

import numpy as np
from scipy.ndimage import uniform_filter

# Define the moving average filter function
def moving_average_filter(image, kernel_size=3):
    """
    Apply a moving average filter to an image.

    Parameters:
    - image: A 2D numpy array representing the image.
    - kernel_size: The size of the kernel (window) for the moving average. Default is 3 .

    Returns:
    - A 2D numpy array representing the filtered image.
    """
    # Apply the uniform filter, which is a moving average filter
    filtered_image = uniform_filter(image, size=kernel_size)
    return filtered_image

# Load the image and convert it to grayscale
image_path = './0183.jpg'
image = Image.open(image_path).convert('L')
image_array = np.array(image)

# Apply the moving average filter with a kernel size of 5
filtered_image_array = moving_average_filter(image_array, kernel_size=5)

# Convert the filtered array back to an image for display
filtered_image = Image.fromarray(filtered_image_array.astype(np.uint8))

# Flatten the original and filtered image arrays for analysis
original_pixel_values = image_array.flatten()
filtered_pixel_values = filtered_image_array.flatten()

# Calculate mean and standard deviation for both original and filtered images
original_mean = np.mean(original_pixel_values)
original_std_dev = np.std(original_pixel_values)
filtered_mean = np.mean(filtered_pixel_values)
filtered_std_dev = np.std(filtered_pixel_values)

# Calculate control limits for both images
original_lower_bound = norm.ppf(0.001, original_mean, original_std_dev)
original_upper_bound = norm.ppf(0.999, original_mean, original_std_dev)
filtered_lower_bound = norm.ppf(0.001, filtered_mean, filtered_std_dev)
filtered_upper_bound = norm.ppf(0.999, filtered_mean, filtered_std_dev)

# Identify outliers for both images
original_outliers = np.where((original_pixel_values < original_lower_bound) | (original_pixel_values > original_upper_bound))
filtered_outliers = np.where((filtered_pixel_values < filtered_lower_bound) | (filtered_pixel_values > filtered_upper_bound))

# Plotting
fig, axs = plt.subplots(2, 2, figsize=(15, 12))

# Original grayscale image and its control chart
axs[0, 0].imshow(image_array, cmap='gray')
axs[0, 0].set_title('Original Grayscale Image')
axs[0, 0].axis('off')
axs[1, 0].plot(original_pixel_values, marker='o', linestyle='', markersize=0.5, color='blue')
axs[1, 0].plot(original_outliers[0], original_pixel_values[original_outliers], marker='o', linestyle='', markersize=0.5, color='red')
axs[1, 0].axhline(y=original_mean, color='r', linestyle='-', label='Mean')
axs[1, 0].axhline(y=original_lower_bound, color='g', linestyle='--', label='Lower Bound')
axs[1, 0].axhline(y=original_upper_bound, color='b', linestyle='--', label='Upper Bound')
axs[1, 0].set_title('Control Chart for Original Image')
axs[1, 0].legend()

# Filtered image and its control chart
axs[0, 1].imshow(filtered_image, cmap='gray')
axs[0, 1].set_title('Moving Average Filtered Image')
axs[0, 1].axis('off')
axs[1, 1].plot(filtered_pixel_values, marker='o', linestyle='', markersize=0.5, color='blue')
axs[1, 1].plot(filtered_outliers[0], filtered_pixel_values[filtered_outliers], marker='o', linestyle='', markersize=0.5, color='red')
axs[1, 1].axhline(y=filtered_mean, color='r', linestyle='-', label='Mean')
axs[1, 1].axhline(y=filtered_lower_bound, color='g', linestyle='--', label='Lower Bound')
axs[1, 1].axhline(y=filtered_upper_bound, color='b', linestyle='--', label='Upper Bound')
axs[1, 1].set_title('Control Chart for Filtered Image')
axs[1, 1].legend()

# Display the plots
plt.tight_layout()
plt.show()

len(original_outliers[0]), len(filtered_outliers[0])

(3482, 3207)